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ABSTRACT 


The representation of unsteady aerodynamic flow fields in terms of global aerody namic modes 
has proven to be a useful method for reducing the size of the aerodynamic model over those 
representations that use local variables at discrete grid points in the flow field. Eigenmodes 
and Proper Orthogonal Decomposition (POD) modes have been used for this purpose with 
good effect. This suggests that system identification models may also be used to represent 
the aerodynamic flow field. Implicit in the use of a systems identification technique is the 
notion that a relative small state space model can be useful in describing a dynamical system. 
The POD model is first used to show that indeed a reduced order model can be obtained 
from a much larger numerical aerodynamical model [the vortex lattice method is used for 
illustrative purposes] and the results from the POD and the system identification methods are 
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then compared. For the example considered, the two methods are shown to give comparable 
results in terms of accuracy and reduced model size. The advantages and limitations of 
each approach are briefly discussed. Both appear promising and complementary in their 
characteristics. 


T ]\T TTJ nm TPTTn\T 
JL ±N .1. ill J U U V-V JL .1. U 

In recent work it has been shown that the use of a modal representation of unsteady aero- 
dynamic flow fields has many advantages [1-4] . Among these is the ability to reduce the 
order of a computational fluid dynamics (CFD) code from many thousands of degrees of 
freedom to several tens of degrees of freedom while retaining essentially the same accuracy 
of representation of fluid forces acting on a wing. Recognizing that the aerodynamic forces 
may be expressed in terms of aerodynamic modes, this suggests that in some cases it may be 
advantageous to use a system identification method to determine a modal representation of 
the aerodynamic forces. For example, this could be done using either numerical data from a 
CFD code or experimental data from a wind tunnel test. Here we explore this possibility and 
compare the results obtained from a system identification procedure and those obtained by a 
direct determination of the aerodynamic modes from a theoretical fluid dynamics model. In 
the present study, a three-dimensional time domain vortex lattice aerodynamic model and 
the Proper Orthogonal Decomposition(POD) method are used to investigate the unsteady 
flow about an oscillating delta wing. The POD modes and the eigenmodes of the vortex 
lattice model are determined starting from a time history of the flow field computed for 
a step change in angle of attack (or plunge velocity). This same set of time histories is 
then provided as the input data for a system identification technique as discussed in [5] to 
determine a dynamic model of the aerodynamic system. The eigenvalues of this identified 
system are then compared to those determined using the POD method and the vortex lattice 
model. Excellent agreement is found between the two results, thus confirming the ability 
of system identification methods to extract useful information from numerically determined 
time histories based upon a theoretical computational model. This also suggests that a sim- 
ilar procedure might be successful that employs numerical data from more complex CFD 
models or from measured experimental data. In the following sections, discussed in order 
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are the Aerodynamic Model, the Proper Orthogonal Decomposition Method, the System 
Identification Method and an Example using numerical data from a simulated time history. 


A t? t> tyvtyt \ lx/TTr 1 A/rnniPT 
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The flow about a cantilevered half-span delta wing is assumed to be incompressible, in viscid 
and irrotational. Here an unsteady vortex lattice method is used to model this flow. Kim [4] 
has also used a similar flow model in his study of the POD or Karhunen-Loeve method for a 
wing of rectangular planform. A typical planar vortex lattice mesh for the three-dimensional 
flow over a delta wing is shown in Figure 1. The plate and wake are divided into a number 
of elements. In the wake and on the wing all the elements are of equal size, Ax, in the 
streamwise direction. For the present calculations, the spanwise element size, Ay, is chosen 
to be equal to Ax. Point vortices are placed on the plate and in the wake at the quarter 
chord of the elements. At the three-quarter chord of each plate element a collocation point, 
is placed for the downwash, the velocity induced by the discrete vortices is set equal to 
the downwash arising from the step angle of attack of the delta wing. Thus one has the 
relationship, 


«•! ' 1 = Ef im Ky rj +1 , i = 1, km (Al) 

where wf* 1 is the downwash at the ith collocation point at time step t + 1, P ; - is the strength 
of the jth vortex. P is normalized by cU and w by U where U is the freestream velocity. Ap- 
is an aerodynamic kernel function. An aerodynamic matrix equation to determine P can be 
expressed as. 


M|{IT ' 1 + |/-1[l ) ? = |T|{u;}' il (A2) 

where L4] and [B\ are aerodynamic coefficient matrices and [T] is a transfer matrix for 
determining the relationship between the global vortex lattice mesh and the local vortex 
lattice mesh on the delta wing itself. Expressions for A, B and T are given in Ref.jl], 


From the fundamental aerodynamic theory, the pressure distribution on the wing can be 
obtained at the jth point in terms of the vortex strengths. The pressure difference between 



the wing plate upper and lower surfaces at x = Xj, y — pj is given by 


^ '{Xj. yj) - pl‘ 
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(r t+1 (^, Vi) + £U(r t+ 1 (xi,y,) - O.o . Vj ))\ 


(A3) 


The local total lift at y = is obtained by integrating the pressure difference along the 
local chord-line. 

. , , r c local , , 

LiocaM = j G Ap- + (x, yj)dx 

where ci oca i = cy-j/l and l is span length and c is the root chord. The total lift is obtained 
by integrating the local total lift along the span. 

= I LlocMdy ('Yl AH 

J ° J=i 

where A y = l/k n and k n is the number of spanwise discrete elements. The total aerodynamic 
lift coefficient, Cl is defined by 


/ <> 
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(A4) 


where A w is the total wing area. In the present example, A w = ~c 2 and Ay = A 


The numerical model is a simple delta wing configuration with a leading edge sweep of 
45 degrees. We use an aerodynamic vortex lattice model with 120 vortex elements on the 
delta wing (km=kn=15) and 525 vortex elements in the wake (kmm=50). The total number 
of vortex elements is N 645. 


SINGULAR VALUE DECOMPOSITION, PROPER OR 
THOGONAL DECOMPOSITION AND BALANCED 
MODES: A GENERIC DISCUSSION 


Let q(j) be the nth flow variable at some spatial point at some time j where n ----- 1, .... N 
and j -= 1, ..... J. Now form the matrix, (), as 



Again note the total number of time steps is J and the total number of flow variables is N. 
For a typical CFD calculation, J might be 1000 and N might be 10000 or more. Hence N 
is much greater than J . 

Now assume a singular value decomposition of Q i.e. 

[Q\ = [mwr ( 2 ) 


where U is a unitary matrix of dimension N x n and V is also a unitary matrix of dimension 
J x n. One may select n and typically n will be less than J. Note that 

[U] T [U] *> [4xn , [V] T [V] * {I\ nxn (3) 

and S is a diagonal matrix of singular values, i.e. 


o l 


02 


( 4 ) 


On 


Now order these singular values such that 


0\ > o 2 > o n (5) 

Form $, the correlation matrix for the POD method. 

m = = {v}imu\ T [u\m r } T - [ v ][ s ] t p ][ vi t w 

Equation (6) implies that V is the eigenvector of the correlation matrix and the corresponding 
eigenvalues are the squares of the singular values. From Equation (2), one may compute 
(assuming that V is normalized so that the magnitude of each eigenvector is unity) 

[Qm - pmivfiv] - mm (?) 

One may also compute U from Equation (7) and further one may compute Q from a 
knowledge of U, V and the singular values using Equation (2). Usually it is easier to compute 
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Q directly from Equation (1). However the representation of Equation (2) may be useful if 
one chooses to decompose Q such that 

[Q] - ([H][E] 1/2 )([E] 1/2 [E] t ) (8) 

With this decomposition the POD modes are said to be "balanced” and these are often put 
forth as an optimum choice for mode selection. 

If there is a truncation in the singular values, i.e. one chooses n to be less than J 
which is much less than N, then Equation (2) may be written in a reduced form. The 
corresponding reduced form for Q approaches the original Q if the neglected singular values 
or POD eigenvalues are sufficiently small compared to those retained. 

Denoting V as the eigenvector matrix for the correlation matrix of dimension J x n, 
noting that Q is a matrix of N x J, and defining, a, as the new unknowns to be determined 
which are the n modal amplitudes of the POD modes, then one may write the original flow 
variables, q, as 

{</ f! } N x 1 = \Q] N <j[V }j . ' 1 (9) 


Substituting this expression into a generic form of the flow equations of motion, i.e. 

{j ( } {Q(q)} i {/*}« (10) 

and pre-multiplying by the transpose of Q gives a reduced order model in terms of the new 
unknowns, a , where the dimension of the vector a is n x 1 with n chosen to be less than J. 

For simplicity, in Equation ( 10) only a single scalar input, u, is shown. The generalization 
to multiple inputs is clear. If Q(q) in Equation (10) is expanded in a Taylor Series about a 
steady flow solution (the time linearized model corresponds to retaining only linear terms in 
q in the Taylor Series), then a particularly simple and attractive form of the reduced order 
model is obtained. 

There is another interesting case to consider which may arise when experimental data 
rather than numerical data from a CFD code are used to construct a reduced order model. 
In this case the number of flow variables that are observed or measured, N, will be relatively 
small and typically N will be less than J, the total number of time steps for which data are 



obtained. Formally the calculation still goes through, but now the number of flow variables 
modeled is much smaller than for a CFD code. Ideally these flow variables would be related 
to the amplitudes of the dominant modes of the flow. 


SYSTEM IDENTIFICATION 

To conform to the usual formalism of system identification theory, we consider a discrete 
time state space model [4]. Note the vortex lattice How model naturally appears in this 
form. For continuous time models, the relationship between discrete and continuous time 
representations must be carefully considered. Formally the procedure is similar, however. 


Assume that the flow field can be described by the discrete-time state-space model, 

x(k i 1) = Ax(k)+ Bu(k) (11) 

y(k) = Cx(k') T Du(k) 

where x(k ) is a state vector of n x 1 (n order of the system), u(k) is an input force vector 
of 7 x 1 ( 7 =number of input), y(k ) is an output measurement vector of f m x 1 (m number 
of output). The matrices A, B. C, and D, are the state matrix, input matrix, output matrix, 
and direct transmission matrix, respectively. 

Let Ui( 0 ) = l(i = 1,2, 7 ) and m(k) = 0 (k = 1,2,....) be substituted into Eq.(ll). 

When the substitution is performed for each input element »,-(()) = 1 (z = .1,2, 7 ), the 

results can be assembled into a pulse response matrix Y with dimension rn by r as follows: 

Y 0 = D, l'i = CB, Y 2 = CAB , , Y k = CA k ~ l B (12) 

The constant matrices in the sequence are known as Markov parameters (See [4]). 

Assume that the scalar quantity q^k) in the matrix Q in Eq.(l) is the pulse response of 
the ith flow variable at the time step k corresponding to a pulse input. A column vector Y 
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can be formed by 


A = C B = 

’ ?i(l) ' 
Q2 ( 1 ) 

, Y 2 = CAB = 

' ?i(2) ’ 

f/2(2) 

, ....Y k = CAA 1 B = 

Qi{k) 

<l2\k) 


Qrn( 1) 


. <bn( 2) . 


Qm ( t ) 


System identification begins by forming the generalized Ilankel matrix II { 0), composed 
of the Markov parameters 

11(0) ■■■ [ V, Y 2 .... Y, , ] (14) 

This matrix is identical to the matrix Q in Eq.(l) with the absence of the last column. Here 
as assumption has been made, i.e., rn > l > n. For the other cases where the number of 
measurement sensors is sufficiently less than the order of the system, i.e., I > n > rn, a 
different Hankel matrix can be formed for system identification (see Chapter 5 in [4]). The 
fundamental rule is that the Hankel matrix must be formed such that its rank is larger than 
the order of the system to be identified. 


In theory, the Hankel matrix H (0) and the state-space model are related by 


H( 0) = [ Y, Y 2 .... ] 

= [ CB CAB .... CA l ~ 2 B ] 

- C[B AB .... A 1 2 /i j (15) 


To determine A, B, C, first decompose the matrix //(()) by using singular value decomposition 
to yield. 


H( 0) 


unv 


\Ur U t] 


y n 

^n—r 


o„ 




[V r V d' 


r, 7 ' i v,y ( \T 

U r E r Vf 

Pr sy 2 ][sy 2 K] 


( 16 ) 


Comparison of Eqs.(15) and (16) established the following equalities 


[U r 


y- 1/2 

/-J r 
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[T} r /2 v r ] • : [ B A A 1 - 2 ] (17) 

The equalities are not unique, but they are balanced because both share the matrix E equally. 
The matrix B is then determined by 

B = The first column of [Ey 2 V)] (18) 

To determine the state matrix /l, another Hankel matrix must be formed, i.e., 

# (1) = [ Y 2 V) .... Yi } (19) 

This matrix is formed by deleting the first column of 11(0) and adding a new column at the 
end of the matrix. As a result, Hi 1) has the following relationship with the system matrices 
A, B, and C 


Substituting Eq.(17) into Eq.(20) thus yields 


//(l) - [ Y 2 Y 3 .... Yi 

- [ CAB C A 2 B .... CA l 1 B 
CA[b AB A 1 2 B ] 


A - \L ! r Y-i "l 1 //( l)[>:y 2 Vv] + 
= yr l/ 2 U T H( 1)F t E _1/2 

r u r \ 1 : r / ~ J r 


( 20 ) 


( 21 ) 


The symbol f means pseudo-inverse. The orthonormal property of U and V shown in Eq.(3) 
has been used to derive Eq.(21). 


Assume that the state matrix A of order r has a complete set of linearly independent 
eigenvectors ipuib, ,Vv with corresponding eigenvalues A ( , A 2 , , X r which are not nec- 

essarily distinct. Define A as the diagonal matrix of eigenvalues and T as the matrix of 


eigenvectors, i.e., 


A 


A 


A 2 


( 22 ) 


A, 
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and 


T 


W\ 



(23) 


The identified A, B, and C can then be transformed to become A, T l B, and (74/. The 
diagonal matrix A contains the information of modal damping rates and damped natural 
frequencies. The matrix 'b - 1 B defines the initial modal amplitudes and the matrix 6Tb the 
mode shape at the sensor points. AH the modal parameters of a dynamic system can thus 
be identified by the three matrices A, and Ctb. The desired modal damping rates and 
damped natural frequencies are simply the real and imaginary parts of the eigenvalues A e , 
after transformation from the discrete-time domain to the continuous-time domain using the 
relation 

A c = ^~Log{ A) (24) 


AN EXAMPLE 

1. POD model: 


In this example J = 100 time steps and Equation (A2) is marched in time to determine F 
for a unit step change in angle of attack. Identifying F with the flow variables, q, we form Q 
and T. See Equations (1) and (6). From Eq.(6). one obtains the POD eigenvector V and the 
corresponding POD eigenvalues of E T E. We normalize the POD eigenvectors V such that 

NfN] - [/] 

The POD eigenvalues calculated are shown in Table 1. Examining these 100 eigenvalues, 
one finds the first few modes are the most important. Also found are some very small 
negative eigenvalues because of the numerical roundoff error in the eigenvalue calculation. 
Note that the POD eigenvalues and eigenvectors are not the eigenvalues and eigenvectors of 
P per se. The latter maybe determined from Eq.(A2) or approximately from the reduced 
order model. 
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( 25 ) 


From Eq.(2), we can determine the unitary matrix U, 

[U] = [QHV1P]- 1 

with a normalization such that [U] 7 [U] = [/]• From Eq.(A2), we have 

{!'}' 1 + I A | ] \n\{\'Y [A\- ] [T}{w} t+l 


From Eq.(9), We have 

{!'} |G]jvxj[V]j „{«}„. 

Substituting Eq.(25) into Eq.(9), we have 


(26) 


{F} = [V 


a f 


(27) 


Substituting Eq.(27) into Eq.(26), we have 


(28) 


Pre-multiplying by the transpose of [L r ][ E| gives a reduced order model in terms of the new 
unknowns, {a}. 


(P]P]) T [r/][E]{a} ,+1 + (lt/](S]f[4]-'[BHe)[S]{a}‘ = ([t/][E]) T [4P[r|M‘ +1 (29) 


Thus a new three-dimensional time domain vortex lattice aerodynamic equation is given 


by 


where 


and 


Wnxl r [^nxrii^lnxl ~~ [■^JnxJV'{'^}jVxl 

[ p \ - (PH ed-HMPDWHWP] 

m = (p] r p])- i ([tf]p]) r [A]- i [Ti 


(30) 


Generally, n << N, and thus a reduced order time domain equation is obtained, for the 
present example. N = 645 and n < 20. Once {a} is determined from Eq.(30), the {F} 4 and 
{r} 4+1 can be determined from Equation (26) and the total lift coefficient of the delta wing 
is calculated using Eq.(A4). 



Figure 2 shows the total lift coefficient of the delta wing vs nondimensional time, r, as 
determined from the full vortex lattice model i.e , using Eq.(A2) with N = 645. 

Figure 3 shows the total lift coefficient of the delta wing vs nondimensional time, r, as 
determined from the reduced order model for various sizes of the reduced order model, i.e, 
n ------- 1. 5 and n ----- 20, using Eq.(30). 

The reduced order model results are very close to the result using the full vortex lattice 
model for the steady state, r — » oo even for n = 1. These results are consistent with those of 
Table 1. As seen from this table, the POD eigenvalue of the first mode is significantly larger 
than the others. So there is a single dominant mode and all others have a relatively small 
contribution to the total lift coefficient. However more POD modes are required to capture 
the dynamics of the flow model for shorter times, say r < 1, as shown in Figure 3. 


The aerodynamic eigenvalues of the vortex lattice model may be determined from Eq.(26) 
and an approximation to these eigenvalues may be determined from Eq.(30) by setting the 
right hand sides to zero and determining the condition for non-trivial solutions in the usual 
way. Results from these calculations will be compared to those obtained from the system 
identification method in the following discussion. 


2. System identification: 

Using the time history for the vortex strengths of the vortex lattice model at each grid point, 
the system identification procedure described earlier was used to construct a dynamical model 
of this system. By using all the information of the original model we can be sure that we will 
reproduce these same time histories to some specified accuracy and that is indeed the case. 
Within plotting error, the results for the time history of the original vortex lattice model 
and that of the identified model are identical. 

More interesting is the number of states required in the system identification model to 
achieve this accuracy and how well the eigenvalues of the identified model agree with those 
of the original model. It turns out that the number of states in the identified model is about 
40, i.e. of the same order as the number of states required in the POD model to reproduce 
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the essential dynamics of the original vortex lattice model. 


It has been shown in previous work that the dominant eigenvalues of the vortex lattice 
model may be represented by the POD modes. In Figure 4 we show a comparison between 
the eigenvalues determined from the full vortex model and those determined from the system 
identification model. As is seen the agreement between the two results is excellent for the 
corresponding eigenvalues. The identified model in fact has represented the eigenmodes in 
the dominant branch of the eigenvalue distribution of the original vortex lattice model and 
done so with a high degree of accuracy. 

The dominant branch of the eigenvalue distribution corresponds to the eigenmodes with 
a spanwise pressure distribution similar to that for the lift on a wing at steady angle of 
attack. Other branches of the eigenvalue distribution correspond to higher spanwise modal 
forms. 


A system identification (SID) model for an unsteady aerodynamic flow has been created 
for several wing motions or gust excitations and corresponding aerodynamic responses. These 
models were derived from numerical simulations using a vortex lattice (VL) model for the 
delta wing with 55 vortex elements on the wing and 300 vortex elements in the wake (e.g. 
km=kn=10, kmm =40 as shown in Fig. 1). In each case, the flow about the wing is excited 
by a certain type of prescribed downwash at the wing grid points, w(t). The numerical VL 
model produced vortex strengths at the grid points of the wing and in the wake, P(£), and 
the corresponding pressure distribution on the wing, p(t). These data were then used as 
input for the SID code. 


The excitations to the flow that have been considered are 

1) step angle of attack, w(t) -- const for t > 0, 

2) sharp edge gust, w(t — x/u) = const for (t — x/u) > 0, 

3) frozen gust of changing frequency, 

wit — x/u) ------ const x sin(a;(£ — x/u) 2 /‘IT) 


such a gust is sometimes called a swept gust, and finally 



4) random downwash (w at each grid point and at every time step is a random number). 


First, the SID model was used to reproduce the time histories of the vortex strength 
(on the wing and in the wake) for 100 time steps. Results obtained from the SID model 
were almost identical to those of the original VL model. The plots of such time histories are 
not provided here, because one would not be able to see the difference in the time series of 


the original data and the time series obtained via SID. However, to quantify the differences 
between VL and SID outputs an error, 5, defined as 5 100%|Q — y|/|Q| was used, where 

the norm |X| is defined to be the largest singular value of X. In this case, both Q (VL) and 
y (SID) are 100 x 355 rectangular matrices (100 time steps and 355 degrees of freedom). For 
the above mentioned aerodynamic problems, and given the vortex strengths everywhere, the 
error was less than 0.05%. Moreover, using the VL time series for the vortex strengths only 
on the wing, the SID model reproduced them with almost the same accuracy - the largest 
(among the aforementioned cases) error was 0.1%. The number of singular values retained in 
the identified model varied from 60 to 90 for the cases when vortex strengths were provided 
only on the wing and from 60 to 160 for the cases when the SID used vortex strengths both 
on the wing and in the wake. (The number of singular values kept was such that the ratio 
of the largest singular value to the smallest one kept was 10 " 12 .) 


Next the ability of such a SID model to capture the eigenmodes of the original VL system 
and especially the dominant eigenmodes was studied. Here again, all the information on the 
wing and wake was used initially, and then the information only on the wing was used to 
obtain the eigenvalues. The former case is shown in Figures 5 and 6 and the latter in Figure 7. 
Here the original VL eigenvalues are shown as circles and the eigenvalues obtained from the 
SID model are indicated by crosses. In the case of full information, the dominant eigenvalue 
branch for the case of a step angle of attack (Fig. 5a) and also a sharp edge gust (Fig. 5b) 
is captured nicely (36 singular values and 57 singular values were used respectively in the 
SID model). For the case of a swept gust (Fig. 5c) and random downwash (Fig. 5d), the 
SID model produced sets of eigenvalues that do not agree well with the original ones from 
the VL code (164 singular values and 99 singular values were used respectively in the SID 
model). However, for the problem of the random downwash, increasing the length of the 
time history that was used as input to the SID model code from 100 steps to 1000 steps lead 
to a set of eigenvalues that is in very good agreement with the dominant branch as shown 
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in Figure 6 (250 singular values were used). 

When the vortex strengths only on the wing were used for some initial selection of 
parameters in the SID model (see Fig. 7), interestingly the dominant branch of eigenvalues 
is recovered not only for the step angle of attack (Fig. 7a) and the sharp edge gust (Fig. 
7b), but also for the case of the swept gust (Fig. 7c) and random downwash (Fig. 7d), 
even though in the last case the dominant branch identified by the SID model is somewhat 
’’off target”. By looking at the Figure 7 one might conclude that, providing vortex strength 
information only on the wing in these cases is sufficient, for finding the distinct dominant 
eigenvalue branch. However, one should not conclude that using vortex strength information 
on the wing only is actually better than using information about vortex strengths both on 
the wing and in the wake. 


By running the SID code for various parameter choices and for a certain number of 
retained singular values, different lengths of time series, using information on the wing only 
vs information everywhere, etc, the authors have found, firstly, that the set of eigenvalues 
from the SID model could either include the '’correct” dominant eigenvalue branch or could 
be a different set of eignvalues (which, by the way, does not prevent the SID model from 
reproducing the original time series with excellent agreement); secondly, it was concluded 
that if there is a distinct branch of eigenvalues found from the SID model, it is most likely 
to be similar to the dominant branch in the original VL model. 


Now another case will be considered. In the context of wing tunnel or flight experiments, 
one is also interested in what can be identified from the measured data such as pressure. 
So using numerically obtained pressure values over the wing, the SID model was used to 
reproduce pressure time histories. Just as in the case of vortex strength data on the wing 
only, 55 pressure time histories of 100 time steps were input to the SID code. Using 30 to 
40 singular values, the error between the original and identified time histories of pressure 
remained less than 1 %. 

One could argue that it may be impractical to use such a large amount of sensors on the 
wing to apply SID models. Thus, next only a portion of the pressure data on the wing was 
used. Different sensor positions were considered. A discussion of the various downwash and 
sensor locations on such a wing is ommitted here. For the current discussion, leading edge 
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elements from the root to the tip of the wing were taken, as they appear to he a very good 
choice for system identification purposes. In Figure 8, results for the pressure time history at 
the tip element of the wing are shown for the case when 10 measurements (still numerically 
simulated, of course) of the down wash and pressure time history (100 time steps) along the 
leading edge were supplied to create the SID model. Using from 50 to 90 singular values for 
the four aerodynamic loading problems, the previously defined error ranged from 0.002% to 
1 , 8 %. 

However, using pressure data only, no success has yet been realized in approximating the 
eigenvalues of the original VL model. This is a subject of continuing investigation. 


CONCLUDING REMARKS 


It has been shown that a modal representation of the aerodynamic forces acting on an 
oscillating wing has several attractive features that, may be exploited. As has been shown 
earlier, starting from a numerical model of the fluid, here taken to be a vortex lattice model, 
the fluid model may be efficiently represented, in terms of Proper Orthogonal Decomposition 
(POD) modes. However it was also shown that using standard system identification tech- 
niques for linear dynamical systems, an identified system model may also be constructed 
using data from the numerical model. Both the POD modes and the system identification 
method lead to substantially reduced order models of the flow field of comparable dimen- 
sions. And therefore both provide highly efficient computational models relative to the 
original [vortex lattice] numerical model. 

In principle the POD method may also be used for nonlinear dynamical systems. How- 
ever much work has yet to be done to implement this capability in the context of unsteady 
aerodynamic models. The system identification methodology for nonlinear dynamical sys- 
tems is still a research frontier and no general methods are yet available. On the other hand 
the systems identification approach may be useful when only partial flow field data are avail- 
able, e.g. from a wind tunnel experiment, and this aspect is currently under investigation 
by the authors. 


Finally it is noted that although a particular numerical model of the fluid was used in the 
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present paper, i.e. a vortex lattice model, no especial difficulty is expected in extending these 
methods to other Computational Fluid Dynamics (CFD) models. Indeed the beneficial use 
of POD for a variety of CFD models has been shown in [1-3]. Since the system identification 
model only requires time histories or Fourier Series (frequency domain) data from such a 
CFD model, the extension to any time linearized CFD model is straightforward. Indeed 
even if the original CFD model is not time linearized per se, the code can be run for small 
wing motions to simulate numerically time linearized conditions. 
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Figure 1: Numerical grid for delta wing-plate using vortex lattice aerodynamic model 



Figure 2: The total lift coefficient of the delta wing vs nondimensional time, r, for all mode 





Figure 3: The total lift coefficient of the delta wing vs non-dimensional time, r, for n = 1, 5 
and n = 20 



0 2 4 6 8 

Frequency, Hz 


Figure 4: Eigenvalue solutions of vortex lattice model of unsteady three-dimensional flow 
about a delta wing. Eigenvalues determined directly from the vortex lattice model: o Dom- 
inant branch of eigenvalues, • All other eigenvalues; Eigenvalues from system identification 
model: x. 
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Figure 5: Original eigenvalues from VL code (o) and eigenvalues from system ID model 
obtained using vortex strengths data on the wing and in the wake (x): (a) step angle of 
attack, (b) sharp edge gust, (c) swept gust, and (d) random case. 
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Figure 6: Original eigenvalues front VL code (o) and eigenvalues from system ID model 
obtained using vortex strengths data on the wing and in the wake (x) for the random case 
when time series length is increased 10 times to 1000 time steps. 
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Figure 7: Original eigenvalues from VL code (o) and eigenvalues from system ID model 
obtained using vortex strengths data on the wing only (x) : (a) step angle of attack, (b) 
sharp edge gust, (c) swept gust, and (cl) random case. 







Figure 8: Time series of the pressure at the tip element. Original VL code data (solid line 
with o’s) and resulted time series from system ID model obtained using pressure on only 10 
leading edge elements as input (dotted line with x’s): (a) step angle of attack, (b) sharp 
edge gust, (c) swept gust, and (d) random case. 
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